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FOREWORD 


This  work  was  sponsored  under  the  auspices  of  the  Naval  Surface  Warfare  Center’s  long-term 
study  program.  This  program  allows  employees  the  opportunity  of  continued  academic  study  for  the 
period  of  1  year.  The  study  was  conducted  during  the  summer  and  fall  of  1987  under  the 
aforementioned  program.  The  purpose  of  this  study  was  to  develop  an  analytical  approach  for  the 
prediction  of  underwater  bubble  collapse  and  its  interaction  with  the  free  surface  or  a  solid  surface. 
The  report  describes  the  solutions  for  axisymmetric  problems.  The  approach  is  shown  to  be  in  good 
agreement  with  available  experimental  data  and  other  analytical  approaches. 
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1.  INTRODUCTION 


An  underwater  bubble  can  exhibit  dynamic  behavior  similar  to  a  mass  spring  system.  For  an 
explosion  bubble,  the  initially  high  internal  pressure  will  drive  the  surrounding  fluid  outward  at  a 
decreasing  rale.  Due  to  the  fluid  inertia,  the  bubble  will  overexpand,  dropping  its  internal  pressure  to 
only  a  small  fraction  of  the  surrounding  fluid’s  hydrostatic  pressure.  Then  the  bubble  will  begin  to 
contract  at  an  increasing  rate.  The  contraction  is  propelled  by  the  fluid's  surrounding  hydrostatic 
pressure.  This  process  will  continue  until  the  increasing  bubble  pressure  abruptly  reverses  the  process. 
Thus,  the  elastic  properties  of  the  gas  and  water  provide  the  conditions  for  an  oscillating  system.  This 
process  is  termed  bubble  pulsation. 

During  the  contraction  phase  of  the  bubble  pulsation,  variations  in  the  surrounding  fluid  pressure 
can  cause  the  bubble  surface  to  contract  at  different  rates.  These  pressure  variations  can  occur  due  to 
bubble  migration,  a  nearby  object,  the  natural  gradient  of  pressure  with  depth  in  a  fluid,  and  other 
disturbances.  These  external  pressure  variations  will  cause  a  higher  rate  of  inward  fluid  acceleration 
near  a  portion  of  the  bubble  surface.  The  portions  of  the  bubble  surface  near  the  higher  fluid 
acceleration  will  form  a  jet  near  the  end  of  the  contraction  phase  of  the  bubble  pulsation.  This  jet  can 
vary  greatly  in  size  and  momentum,  depending  on  the  factors  leading  to  its  formation.  This 
phenomenon  is  termed  bubble  jetting.  The  primary  purpose  of  this  study  was  the  development  of 
analytical  techniques  that  could  be  used  to  study  bubble  pulsation  and  bubble  jetting. 

Much  of  the  interest  in  underwater  bubble  dynamics  has  focused  on  the  behavior  of  cavitation 
bubbles.  Cavitation  bubbles  have  applications  affecting  military  strategy,  propeller  and  turbine  blade 
damage,  and  underwater  acoustics.  Collapsing  cavitation  bubbles  can  generate  tiny  regions  where 
temperatures  arc  in  the  thousands  of  degrees  with  tremendous  pressures  and  f!  d  velocities.  A 
number  of  authors  (Gucrri,  Lucca,  and  Prospcrctti  1982;  Lezzi  and  Prospcretli  1987;  Prosperctti  1982a, 
1982b,  1984.  1986,  1987;  Prosperctti,  Crum,  and  Commander  1988;  Prospcrctti  and  Jones  1984; 
Prospcrctti  and  Lezzi  1986;  Chahinc  1982,  1977;  Chahinc  and  Bovis  1983;  Chahinc  and  Gcnoux  1983; 
Chahinc,  Gcnoux,  and  Liu  1984;  Chahinc  and  Sivian  1985;  Gcnoux  and  Chahinc  1984;  Johnson  ct  al.; 
Taib  1985;  Chapman  and  Plcssct  1970,  1972)  have  demonstrated  that,  in  some  instances,  bubble 
behavior  can  be  accurately  predicted  with  good  accuracy  utilizing  numerical  methods. 


> 

The  focus  of  this  numeric  study  is  the  dynamic  expansion  and  collapse  of  large  bubbles,  such  as 
those  created  by  an  underwater  explosion  in  the  free  field,  in  the  vicinity  of  a  free  surface,  and  near  a 
solid  body.  Bubble  collapse  is  influenced  by  a  number  of  physical  factors.  The  factors  include  the  * 

gradient  of  hydrostatic  pressure  as  a  function  of  depth,  bubble  migration,  and  a  bubble’s  proximity  to 
nearby  solid  bodies. 

The  effect  of  the  pressure  gradient  on  jetting  is  far  more  influential  for  a  large  bubble  than  a 
cavitation  bubble.  In  the  former  case,  the  hydrostatic  pressure  at  the  bottom  of  the  bubble  surface  can 
be  twice  as  high  as  the  hydrostatic  pressure  at  the  top  of  the  bubble.  Contrarily,  bubble  migration 
leads  to  bubble  jetting  in  both  large  and  small  bubbles.  In  general,  bubbles  are  driven  upward  by  their 
own  buoyancy  force.  Consequently,  the  higher  velocity  of  the  flow  at  the  head  of  the  bubble  leads  to 
a  reduction  in  hydrostatic  pressure,  while  the  stagnant  flow  along  the  lower  surface  of  the  bubble  leads 
to  an  increase  in  the  hydrostatic  pressure.  This  phenomenon,  accompanied  by  the  dynamic  collapse  of 
the  bubble,  will  typically  yield  a  bubble  jet  near  the  end  of  the  bubble  collapse  phase.  Further,  the  t 

proximity  of  the  bubble  to  a  solid  boundary  can  also  greatly  influence  the  bubble  jetting  behavior 
(Figure  1 ).  When  a  bubble  pulsates  near  a  body,  the  fluid  near  the  object  must  move  faster  than  the 
fluid  around  other  portions  of  the  bubble.  The  result  is  a  decrease  in  the  fluid's  pressure  near  the 

boundary.  This  pressure  drop  in  the  fluid  near  the  body  will  generate  bubble  migration  toward  the  1 

body.  When  the  bubble  comes  in  contact  with  the  body  during  the  expansion  phase  of  its  period,  the 

bubble  can  have  a  very  high  water  velocity.  This  is  in  part  due  to  the  dynamic  nature  of  the  bubble 

growth.  As  the  bubble  expands,  the  internal  pressure  can  drop  to  only  a  small  fraction  of  the  normal 

hydrostatic  pressure  in  the  surrounding  fluid  field.  The  result  is  the  jetting  of  the  bubble's  lower 

surface  into  the  body.  Finally,  the  combination  of  the  various  effects  influencing  bubble  collapse  and 

jetting  can  act  in  conjunction  with  one  another  to  increase  jetting  or  oppose  one  another  to  reduce 

jetting.  The  present  method  described  in  this  report  allows  the  opportunity  for  these  behaviors  to  be  i 

studied  separately  or  in  combination  with  one  another. 

The  present  study  utilizes  an  axisymmctric  boundary  integral  formulation  which  has  application  for 
a  variety  of  geometric  configurations.  The  method  is  computationally  efficient  due  to  its  dependence  * 

on  only  the  boundaries'  value  for  the  velocity,  the  potential,  and  the  potential's  first  derivative.  This 
approach  avoids  the  necessity  of  solving  Laplace’s  equation  in  the  entire  domain  occupied  by  the 

> 


2 


liquid  which  would  be  computationally  intensive.  The  resulting  range  in  computation  time  can  vary 
from  minutes  for  free  field  calculations,  to  less  than  one-half  hour  for  computations  involving  free 
surfaces  or  a  boundary.  The  computational  device  serving  as  a  basis  for  these  comparisons  is  the 
VAX  8800.  To  illustrate  the  performance  and  accuracy  of  the  method,  a  number  of  comparisons  were 
made.  These  comparisons  included:  exact  vs.  predicted  solutions  for  spherical  bubbles;  experimental 
bubble  radii  and  period  comparisons;  and  bubble  jet  tip  velocities  vs.  compressible  flow  theory 
calculations  using  the  hydrocode  PISCES.  In  these  cases,  the  method  is  shown  to  be  in  good 
agreement  with  available  data.  Further,  a  number  of  calculations  very  near  the  free  surface  are 
performed  showing  the  idealized  formation  of  the  water  spray  dome.  Using  this  free  surface 
calculation,  the  limitations  of  the  methodology  are  discussed. 

The  growth  and  collapse  of  large  bubbles,  such  as  those  created  by  underwater  explosions  near 
nonaxisymmetric  underwater  structures,  are  also  of  interest.  In  this  scenario,  the  combination  of 
gravity-induced  collapse  and  collapse  caused  by  migration  and  proximity  to  nearby  structures  can 
result  in  extremely  high  velocities  of  the  bubble  reentry  jet  tip.  Analysis  of  bubble-structure 
interaction  is  generally  three-dimensional  in  nature  due  to  the  effects  of  gravity  and  the  orientation  of 
the  structure.  Solution  of  this  problem  is  briefly  discussed  in  Section  2. 

2.  MATHEMATICAL  FORMULATION 

In  this  mathematical  formulation  model,  the  flow  induced  by  one  or  more  bubbles  is  considered. 
The  fluid  occupies  a  domain  ft  bounded  by  the  bubble  surfaces,  S^,  solid  boundaries.  Sr,  and  a  surface 
at  infinity,  Sx  (Figure  2).  The  fluid,  contained  in  domain  ft,  is  considered  incompressible  and 
inviscid. 

In  the  case  of  underwater  explosive  detonations,  a  chemical  reaction  converts  the  original  material 
into  a  gas  at  very  high  temperatures  and  pressures.  The  temperatures  can  be  of  order  3.000°  C  with 
pressures  of  order  50,000  atm  (Landau  and  Lifshitz  1959).  The  result  of  the  initial  detonation  is  a 
compression  or  shock  wave  being  emitted  into  the  fluid  field,  followed  by  the  dynamic  expansion  and 
contraction  of  the  gas  bubble.  After  the  release  of  the  shock  wave,  which  is  an  early  time 
phenomenon,  the  speed  of  the  bubble  surface  remains  an  order  of  magnitude  smaller  than  the  speed  of 
sound  in  water.  Therefore,  imposing  an  incompressibility  condition  on  the  fluid  is  deemed  valid  for 
this  flow  condition. 


Using  the  bubble  radius  as  a  characteristic  length  in  the  expansion  or  contraction  phase  of  the 
bubble  pulse,  a  Reynolds  number  can  be  calculated.  The  Reynolds  number  for  an  underwater 
explosion  is  found  to  be  high  through  most  of  the  bubble  growth  and  collapse.  Recalling  that  the 
Reynolds  number  is  a  ratio  of  inertial  to  viscous  forces,  it  is  clear  that  neglecting  viscid  terms  in  the 
conservation  of  momentum  equation  detracts  very  little  from  the  problem  solution.  Furthermore,  under 
the  assumption  that  the  flow  is  inviscid  and  irrotational  flow,  the  velocity  field  can  be  found  from  the 
gradient  of  a  potential.  In  the  case  of  a  cavitation  bubble,  similar  assumptions  were  shown  to  be  valid 
by  Guerri,  Lucca  and  Prosperetti  (1982). 

The  mathematical  problem  described  above  can  be  stated  in  the  following  form.  The  divergence 
of  the  velocity  vector  is  zero  for  an  incompressible  fluid.  This  expression  is  found  from  the 
conservation  of  mass  equation  and  can  be  written  as 

V  •  u  =  0.  (1) 

For  irrotational  flow,  the  curl  of  the  velocity  vector  is  zero  and  can  be  written  as 

V  x  u  =  0.  (2) 

This  expression  is  also  satisfied  when  the  velocity  u  equals  the  gradient  of  a  potential  which  is 

u  -  V<p .  (3) 

Combining  Equations  1  and  3  yields  Laplace’s  equation.  Laplace’s  equation  is  satisfied  in  the  domain 
£i  occupied  by  the  fluid  and  can  be  written  in  terms  of  the  potential  as 

V2<t>  =0.  (4) 

The  boundary  conditions  necessary  for  the  solution  of  Laplace’s  equation  in  the  domain  £2  can  be 
stated  as 


dx 

dt 


V  0  on  Sh 


(5) 
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=  V<J>  ■  n  =  0  on  Sl, 
dn 

where  x  is  the  distance  to  the  boundary  from  a  consistent  point  of  reference. 

The  conservation  of  momentum,  or  Navier  Stokes  equation  under  the  assumption  of 
incompressible  fluid  and  irrotational  flow,  can  be  written  in  the  form 

Du  _  -Vp 
~Dt  ~  ~ p~  +8' 


(6) 


(7) 


(8) 
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where  Du  represents  the  total  derivative,  P  the  pressure,  p  the  density,  and  g  the  gravity.  Using  the 
identity: 


Vu2 

u  -  Vu  =  _r_  - «x  V  x  u 


(9) 


I 


and  the  relations  in  Equations  2  and  3,  the  equation  for  conservation  of  momentum  reduces  to  the 
form  given  by  Euler  (Landau  and  Lifshitz  1959)  as 


3$ 

~SF 


(10) 


with  PR  representing  the  pressure  on  the  surface  of  the  bubble  and  P^  the  pressure  at  infinity.  The 
pressure  can  be  determined  by  the  choice  of  a  reference  point.  For  example,  if  the  point  of 
reference  at  infinity  is  the  free  surface  of  the  water,  P^  =  1  atm  and  z  is  the  depth  of  the  bubble 
segment.  PR  can  be  found  from  the  adiabatic  pressure  balance  equation  given  by  Cole  (1948)  as 
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(11) 


where  W  is  the  weight  of  the  explosive  in  grams  of  TNT  equivalent,  V  is  the  volume  in  centimeters 
cubed,  and  the  pressure,  PR,  is  given  in  kilobars.  The  o  term  is  added  to  account  for  the  effects  of 
surface  tension.  For  a  cavitation  bubble,  a  similar  adiabatic  expression  can  be  used. 

The  formulation  strategy  for  the  time  integration  of  the  equation  of  motion  is  straightforward.  If 
the  bubble  surface  Sb,  internal  pressure,  and  the  value  of  the  potential  and  its  first  derivative  on  the 
bubble  surface  are  known,  Equations  7  and  10  can  be  solved  while  marching  through  time.  In  order 
to  solve  for  the  potential’s  first  derivative  on  the  bubble  surface,  Laplace’s  equation  (Equation  4)  must 
be  solved.  This  can  be  accomplished  using  the  boundary  conditions  given  in  Equations  5,  6,  and  7 
and  using  a  boundary  integral  method.  The  starting  point  of  this  solution  to  Laplace’s  equation  makes 
use  of  Green’s  second  identity  which  is 

jv  [<t>  V2  v  -  =  J5[<t>V\|f  •  n-\|»V<|>  ■  n]dS  (12) 

with  y  being  a  function  of  choice.  The  V  represents  integration  over  the  volume  and  5  represents 
integration  over  the  surface.  The  most  obvious  choice  for  \y  is 


V  =  - 

I*  *  *1 


(13) 


which  is  the  solution  of  the  Poisson’s  equation  having  the  form 

V2t| /  =  -4rc6(x  -x)  (14) 

where  x  and  x  are  the  reference  and  current  coordinates  of  the  bubble  surface  (Figure  3).  The  choice 
of  Poisson’s  equation  and  the  Dirac  function  are  ideal  because  the  remaining  volume  integral  given  on 
the  left-hand  side  of  Equation  12  reduces  to  a  constant  along  each  bubble  segment. 


6 


► 


* 


» 


a 


» 


> 


» 


> 


» 


( 


> 


> 


* 


With  the  substitution  of  Equations  4,  13,  and  14  into  Equation  12,  the  formulation  reduces  to  a 
surface  integral  with  the  following  form: 


2tc<J>  (x) 


•  n  -  \yV<p  •  n]dS. 


(15) 


» 


I 


The  factor  of  1/2  in  the  left-hand  side  of  the  integration  in  Equation  15  arises  from  the  fact  that  the 
integration  on  the  bubble  boundary  can  only  be  carried  out  around  half  the  point  of  reference  which  is 
the  portion  in  the  domain  fi. 

I 

Equation  15  applies  to  both  three-dimensional  developments  and  the  present  axisymmetric 
formulation.  The  three-dimensional  approach  can  be  accomplished  using  first-order  and  second-order 
quadrilateral  or  triangular  finite  element  techniques.  The  resulting  integration  would  proceed  in  much 
the  same  way  as  in  the  axisymmetric  case.  However,  the  biggest  difference  in  the  three-dimensional  » 

solution  is  in  the  computation  times  required  to  solve  problems.  The  computation  time  for  the  three- 
dimensional  case  would  undoubtedly  be  an  order  of  magnitude  higher  than  the  axisymmetric 
formulation.  Therefore,  the  present  axisymmetric  formulation  affords  some  basic  advantages  in  the 
study  of  bubble  jetting  phenomena.  The  real  drawback  is  the  inability  of  the  present  method  to  tackle  * 

more  complex  interactions  of  bubbles  and  various  underwater  objects  such  as  cylindrical  and  spherical 
bodies  in  nonaxisymmctric  orientations. 

The  axisymmetric  formulation  begins  by  a  suitable  choice  of  coordinate  systems.  The  choice  of 
the  cylindrical  coordinate  system  for  this  case  will  reduce  the  integration  around  the  bubble  surface  to 
an  elliptic  integral.  Using  cylindrical  coordinates,  we  can  write  the  following: 

i 

x  =  ( rcosG ,  rsinQ,  Q  (16) 


for  the  reference  coordinate  system  and 


> 


x  =  (/?,0,Z) 


(17) 


for  our  current  coordinate  system.  Substitution  of  these  cylindrical  coordinates  into 
Equation  15  and  the  rearrangement  of  terms  reduced  the  integration  to  the  following  form: 
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J_  f  r  d<f  p* 
2jtJr  -gjJb 


de 

\-X 


£ 


in  de 


X -X\ 


dr 


(18) 


where  dT  represents  the  integral  on  a  surface  element  (Figure  3)  and  0  the  angle  around  the  bubble 
circumference.  The  integral  about  d6  can  be  reduced  to  an  elliptic  integral  of  the  first  kind  with  a  fair 
amount  of  algebraic  manipulations.  The  elliptical  integral  introduces  a  logarithmic  term  which  is 
singular  for  the  cases  when  the  current  coordinate  and  the  reference  coordinate  coincide.  The  resulting 
singularity  in  the  integration  is  bounded  and  occurs  due  to  the  use  of  the  Dirac  delta  function.  An 
accurate  calculation  of  this  singularity  will  result  in  a  more  stable  and  accurate  bubble  model.  After 
experimentation  with  various  expansions,  a  Gaussian  quadrature  with  weight  log  (!;)  was  found  to  give 
the  best  results  (Wilkerson  1987).  The  quadrature  allowed  the  estimation  of  error  within  acceptable 
limits.  The  final  form  of  the  integration  can  be  written  as 

J-  C  G(x.x)dr  =  <t»(x)  *  -L< t>(x)  f  H(x.x)dr  (19) 

2k  dn  Jr  2*  ir  — 

where  G(x,  x)  and  H(x,  x)  are  functions  of  the  reference  and  current  coordinates  and  an  elliptic 
integral.  Using  a  straight  line  approximation  for  the  axisymmetric  bubble  boundary  (Figure  3),  the 
gradient  of  the  potential  can  be  found  using  the  following  matrix  expression  by  a  Gaussian  elimination 
routine.  The  basic  expression  is 


n 

Y.  aij®j =  b>  ‘ =  ’•2'3-n 

7*1 


(20) 


where 


(21) 


(22) 
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and 


» 


d<b 

h>  =  __  =  const 
1  dn 


(23) 

» 


The  solution  for  d<p/dn  can  now  be  applied  to  the  conservation  of  momentum  equation  and 
marched  through  time.  The  updating  of  line  segments  and  the  momentum  equation  given  in  Equation 
7  is  accomplished  with  a  finite  difference  time  integration  scheme  and  linear  interpolation  between  the 
nodes.  The  first-order  form  of  the  finite  difference  formula  is 


r  i 

+  S* 

2 

L  p  J 

» 


(24) 


However,  higher  order  formulations,  depending  on  multiple  values  of  the  potential,  are  used  after  the 
initial  time  step.  A  similar  expression  is  used  to  update  the  geometry.  These  higher  order 
formulations  do  not  increase  computation  time  significantly.  The  majority  of  the  computation  time 
needed  in  the  calculation  is  spent  doing  the  integrations  in  Equation  19  and  assembling  and  solving  the 
matrix  formulation  given  in  Equation  20.  Therefore,  the  contributing  factors  influencing  computation 
time  are  the  number  of  segments  and  the  number  of  time  steps.  The  amount  of  line  segments 
determines  the  total  number  of  integrations  in  a  time  step  and  also  the  size  of  the  matrix  (Equation  20) 
that  must  be  solved.  For  the  free  field  cases,  32  line  segments  were  found  to  be  sufficient.  However, 
in  the  cases  involving  the  free  surface  or  a  solid  surface,  84  line  segments  were  used.  No  attempt  to 
optimize  the  method  was  made  in  the  present  study.  Contrarily,  the  effect  of  time  step  size  was 
investigated.  In  the  free  field  cases,  any  time  step  less  than  l/450lh  of  the  total  bubble  period  was 
found  to  cause  problems  at  the  beginning  and  at  the  end  of  the  bubble  period.  This  is  probably  due  to 
the  relative  speed  of  the  individual  line  segments  during  these  phases  of  the  bubble  oscillation  period. 
A  formulation  involving  a  variable-sized  time  step  dependent  on  the  highest  speed  of  a  segment  would 
probably  result  in  greatly  reduced  computation  times.  For  the  present,  the  author’s  primary  concern 
was  with  the  validation  and  accuracy  of  the  methodology  and  with  exploration  of  the  limitations  of 
this  formulation.  Consequently,  a  1/1 ,000th  of  the  bubble  total  period  was  used  in  all  cases  examined. 
The  bubble  period  was  estimated  using  standard  empirical  rules  (Swift  and  Dccius  1950)  for 
underwater  explosive  detonations. 
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3.  RESULTS 


The  validation  of  this  methodology  was  initiated  by  attempting  to  solve  simplistic  problems.  The 
problems  were  gradually  changed  to  include  more  complex  phenomena  such  as  bubble  interactions 
with  the  free  surface  and  solid  bodies.  The  problems  and  limitations  are  discussed  with  possible 
solutions  extending  the  method.  What  follows  is  a  discussion  of  various  problem  types,  the  limitations 
of  the  current  approach,  and  possible  advantages  in  the  present  methodology. 

The  starting  point  of  the  validation  was  a  comparison  of  the  boundary  integral  method  to  the  exact 
solution  for  a  spherically  symmetric  bubble  model.  The  conservation  of  momentum  equation  given  in 
Equation  10  was  modified  to  exclude  gravity.  With  gravity  removed  and  in  the  absence  of  other 
bodies,  the  solution  technique  should  duplicate  the  exact  formulation  for  a  spherically  symmetric 
model.  The  conservation  of  momentum  equation  under  spherical  symmetry  and  the  conservation  of 
mass  equation  can  be  integrated  to  yield 


R0 

3 

R2  * 

Y 

Pr  - 

i  - 

R 

o 

3 

L  P  J 

R 

(25) 


where  R  represents  the  current  velocity  of  the  bubble  surface;  R  the  current  radius;  and  PR,  and  p 
(as  given  in  Section  2)  are  the  internal  bubble  pressure,  reference  pressure,  and  fluid  density.  The 
subscript  o  represents  the  referenced  values.  Using  the  boundary  integral  method  for  R,  RQ,  and  Rq, 
the  estimated  value  for  R  was  compared  to  the  exact  solution  given  by  Equation  25.  In  both  bubble 
expansion  and  bubble  contraction,  the  boundary  integral  method  approached  the  exact  solution,  or  the 
limit  of  the  methods  accuracy.  The  accuracy  of  the  method  appeared  to  improve  with  time.  This  is 
caused  by  the  potential’s  initial  condition  being  assumed  zero.  Therefore,  the  method  requires  a 
number  of  steps  to  calculate  a  reasonable  value  for  the  potential.  Furthermore,  this  assumption  is 
particularly  good  in  the  case  of  an  underwater  explosion  where  the  internal  pressure  is  large  by 
comparison  to  the  external  pressure.  The  result  is  that  the  early  bubble  motion  is  dominated  by  the 
large  pressure  differential,  thus  allowing  the  method  to  retain  good  overall  accuracy  under  the  zero 
initial  potential  assumption.  This  accuracy  is  shown  for  both  a  first-order  finite  difference  integrator 
and  a  second-order  time  integrator  in  Tables  1  and  2,  respectively. 


The  convergence  seen  in  Tables  1  and  2  occurs  for  two  reasons.  First,  the  assumption  allowing 
the  potential  to  be  zero  initially  will  require  a  few  time  steps  to  achieve  a  reasonable  level  of 
equilibrium.  This  error  will  decrease  rapidly  at  first  then  gradually  thereafter.  The  second  reason  for 
continued  convergence  is  the  constant  reduction  in  pressure  difference  as  the  bubble  volume  grows. 
This  reduction  in  pressure  difference  results  in  a  decrease  in  the  radial  velocity  of  the  bubble.  Slower 
velocities  under  a  constant  time  step  ate  attributed  to  increased  convergence.  To  ettimate  the  effects 
of  bubble  radial  velocity  on  convergence,  the  method  was  used  with  a  constant  time  step  under 
constant  pressure.  The  constant  pressure  will  result  in  a  constantly  increasing  radial  velocity.  This 
should  allow  an  estimation  of  the  error  attributed  in  time  step  size  vs.  bubble  radial  velocity.  The 
results  are  summarized  in  Tables  3  and  4. 

The  results  showed  that,  for  a  constant  time  step  of  .01  and  velocities  sufficiently  low,  the 
convergence  and  accuracy  were  quite  reasonable.  However,  if  the  bubble  were  allowed  to  iterate  long 
enough,  or  if  the  pressure  differential  was  increased,  the  convergence  reached  a  minimum  and  then 
began  to  diverge  in  accordance  with  bubble  radial  velocity.  By  comparing  Tables  1  and  2  to  Tables  3 
and  4,  the  error  effects  of  a  particular  time  step  at  a  given  velocity  can  be  seen.  In  Tables  1  and  2.  the 
time  step  used  was  .001  for  velocities  of  order  102,  while  Tables  3  and  4  have  a  time  step  of  .01  for 
speeds  of  order  100.  As  expected,  Tables  3  and  4,  which  have  a  smaller  velocity  vs.  time  step  ratio, 
converge  faster  to  a  reasonable  error.  In  as  much  as  the  solutions  converged  to  reasonable  values  in 
both  cases,  it  indicated  that  the  initial  conditions  resulted  in  a  reasonable  cumulative  error.  The  issue 
of  time  step  size  seems  to  be  the  most  crucial  aspect  in  keeping  the  cumulative  error  to  a  minimum. 
This  became  particularly  important  for  bubble  jetting  where  the  velocities  were  as  high  as  1 ,000  ft/s 
for  the  bubble  segments  near  the  nose  of  the  reentrant  jet  tip.  These  high  velocities  lead  to  numerical 
instabilities  which  would  slop  the  analysis.  Therefore,  the  bubble  maximum  velocities  were  compared 
to  the  bubble  period  so  that  a  rule  of  thumb  could  be  developed  for  the  time  step.  It  was  found  that  a 
time  step  of  1/1 ,000th  of  the  bubble  period  would  give  reasonable  results  in  all  cases  investigated.  A 
further  analysis  of  the  time  step  could  lead  to  a  variable  step  size  dependent  on  only  the  maximum 
speed  of  the  fastest  bubble  segment.  However,  no  attempt  to  find  a  suitable  ratio  of  time  step  to 
segment  velocity  has  been  made. 

The  second  phase  of  the  validation  was  to  compare  bubble  periods  and  maximum  bubble  radius  to 
experimental  results.  Swift  and  Dccius  (1950)  offered  data  for  an  abundant  variety  of  explosive  sizes 
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Table  1.  Bubble  Exact  vs.  Calculated  Variable  Pressure 
(First-Order  Time  Integration) 


Radial  Velocity 

Time  Step 

Exact 

Calculated 

Error  (%) 

10 

51.24 

53.10 

3.6 

20 

33.68 

33.77 

0.25 

30 

25.42 

25.37 

0.19 

40 

8.1 

8.18 

— 

9 

Table  2.  Bubble  Exact  vs.  Calculated  Variable  Pressure 
(Second-Order  Time  Integration) 


Time  Step 

Radial  Velocity 

Error  (%) 

Exact 

Calculated 

10 

58.04 

56.67 

2.4 

20 

34.93 

34.94 

0.02 

30 

25.97 

25.97 

— 

Table  3.  Bubble  Exact  vs.  Calculated  Constant  Pressure 
(First-Order  Time  Integration) 


Time  Step 

Radial  Velocity 

Error  (%) 

Exact 

Calculated 

10 

0.09091 

0.09036 

0.60 

20 

0.1944 

0.1943 

0.046 

30 

0.3067 

0.3066 

0.012 
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Table  4.  Bubble  Exact  vs.  Calculated  Constant  Pressure 
(Second-Order  Time  Integration) 


Radial  Velocity 

Time  Step 

Exact 

Calculated 

Error  (%) 

10 

0.0999 

0.1004 

0.48 

20 

0.2048 

0.2047 

0.06 

30 

0.3177 

0.3177 

— 

and  depths  to  be  attempted.  A  number  of  the  shots  presented  in  the  Swift  report  were  selected  at 
random  and  compared  to  the  solution  technique  presented  in  this  report.  The  results  are  presented  in 
Tables  5  and  6.  The  results  generated  in  Table  5  are  of  greater  importance  because  the  methods  used 
to  record  bubble  period  are  more  accurate  than  maximum  bubble  radius  data.  The  accuracy  stems 
from  the  use  of  piezoelectric  gauges  which  can  record  the  pressure  waves  emitted  initially  and  with 
each  successive  bubble  pulse.  The  initial  shock  wave  indicated  the  beginning  of  the  period  and  the 
first  bubble  pulse,  which  occurs  just  after  bubble  minimum,  indicates  the  end  of  the  period.  Due  to 
the  high  speed  of  sound  in  water,  the  method  gives  very  accurate  results.  For  the  data  compared  in 
Table  5,  the  analytical  method  gives  reasonable  results. 

On  the  contrary,  the  maximum  bubble  radii  data  for  the  Swift  report  was  not  as  accurate.  These 
data  were  obtained  from  measurements  taken  from  simultaneous  photographic  records.  As  Swift  and 
Dccius  (1950)  pointed  out,  the  relative  narrow  angle  of  view  from  the  high-speed  cameras  and  the 
problems  with  distortion  and  lighting  made  it  difficult  to  make  bubble  radii  measurements  better  than 
2-4%.  Further,  depth  measurements  were  obtained  from  a  weighted  line  which  could  vary  slightly 
due  to  ocean  currents.  These  factors  are  attributed  to  a  loss  of  accuracy  in  the  bubble  radii 
comparisons  shown  in  Table  6. 

The  accuracy  of  the  method  was  summarized  by  making  a  number  of  radius  vs.  time  plots  in 
comparing  computational  values  to  those  found  in  the  Swift  and  Dccius  (1950)  report.  These  plots  are 
given  in  Figures  4-9.  The  plots  follow  the  growth  and  decay  profiles  of  the  bubble  first  period  quite 
well.  These  plots  represent  results  generated  by  tests  involving  a  number  of  charge  weights  at  varying 
depths. 
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Table  5.  Bubble  Period  (Calculated  vs.  Experiment) 


Shot 

No. 

Charge 

Weight 

(lb) 

Charge 

Dept 

(ft) 

Period 

(ms) 

Calculated 

(ms) 

Error 

(%) 

GIF 

0.651 

343. 

27.40 

26.40 

1.9 

G2F 

0.660 

304. 

29.80 

29.76 

0.13 

G5F 

0.662 

305. 

29.89 

29.70 

0.64 

G6F 

0.658 

304. 

29.92 

29.73 

0.64 

G7F 

0.669 

298. 

30.23 

30.39 

0.50 

G8F 

0.663 

304. 

29.82 

29.80 

0.0 

G9F 

0.651 

304. 

29.64 

29.62 

0.0 

G17F 

0.660 

305. 

29.62 

29.68 

0.20 

G18F 

0.660 

302. 

29.61 

29.92 

1.0 

G20F 

0.660 

538. 

19.64 

18.90 

3.7 

G21F 

0.651 

539. 

19.10 

17.50 

9.0 

G23F 

0.658 

567. 

18.14 

18.30 

0.88 

G70F 

0.660 

503. 

19.90 

20.00 

10.50 

G71F 

0.658 

463. 

21.00 

21.40 

1.90 

G72F 

0.660 

586. 

17.85 

17.80 

0.30 

G73F 

0.655 

576. 

18.16 

18.00 

0.90 

G74F 

0.660 

556. 

18.10 

18.60 

2.80 

G76F 

0.660 

587. 

17.10 

17.38 

0.30 

AVERAGE  1.0% 
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Table  6.  Bubble  Period  (Calculated  vs.  Experiment) 


Shot 

No. 

Charge 

Weight 

(lb) 

Charge 

Dept 

(ft) 

Period 

(ms) 

Calculated 

(ms) 

GIF 

0.651 

343. 

— 

17.0 

G2F 

0.660 

304. 

18.80 

18.49 

G5F 

0.662 

305. 

18.80 

18.48 

G6F 

0.658 

304. 

18.80 

18.47 

G7F 

0.669 

298. 

18.90 

18.70 

G8F 

0.663 

304. 

18.90 

18.52 

G9F 

0.651 

304. 

18.70 

18.40 

G17F 

0.660 

305. 

19.40 

18.46 

G18F 

0.660 

302. 

19.40 

18.53 

G20F 

0.660 

538. 

15.90 

15.00 

G21F 

0.651 

539. 

15.70 

14.40 

G23F 

0.658 

567. 

15.61 

14.70 

G70F 

0.660 

503. 

— 

15.35 

G71F 

0.658 

463. 

16.60 

15.80 

G72F 

0.660 

586. 

15.40 

14.50 

G73F 

0.655 

576. 

15.60 

14.60 

G74F 

0.660 

556. 

15.70 

14.80 

G76F 

0.660 

587. 

15.10 

14.50 

AVERAGE  4.0% 


The  deviations  of  values  obtained  from  theory  found  in  these  figures  is  approximately  of  the  same 
magnitude  as  the  errors  given  in  Tables  5  and  6.  For  deeper  shots,  a  larger  error  occurs  with  the 
radius  calculations.  This  is  believed  to  be  caused  by  a  number  of  factors.  As  pointed  out  earlier,  the 
experimental  data  is  not  as  accurate  for  radius  calculations  as  it  is  for  the  bubble  period.  Further,  the 
initial  conditions  seem  to  play  an  important  role  in  maximum  bubble  radius  calculations.  When  the 
initial  radius  was  decreased,  the  radius  estimates  improved.  This  is  due  to  the  increase  in  hydrostatic 
pressure  at  increased  depths  and  the  effect  it  plays  on  the  overall  pressure  difference.  Additionally,  the 
effects  of  bubble  radius  calculations  are  less  forgiving  for  smaller  charges  than  for  larger  ones. 
However,  in  the  interest  of  consistency,  all  results  calculated  in  this  report  represent  the  same  initial 
conditions  for  internal  bubble  pressure  and  the  potential  on  the  surface. 

All  of  the  data  compared  thus  far  represents  relatively  small  charges.  Therefore,  bubble  periods 
for  a  300-lb  charge  detonated  at  shallow  depths  are  compared  to  results  presented  in  Cole  (1948). 
Figure  10  compares  the  observed  periods  in  Cole  vs.  those  predicted  by  theory.  These  calculations 
include  the  effects  of  a  free  surface  in  the  flow  field.  For  the  shallower  depths.  Figures  11-14  show 
the  growth  and  collapse  of  the  bubble  surface  and  the  effects  of  the  free  surface. 

When  the  bubble  is  closer  to  the  free  surface,  the  effects  of  the  free  surface  become  more  evident. 
This  can  be  seen  in  Figures  11-13  in  terms  of  bubble  migration.  The  migration  is  retarded  due  to  its 
close  proximity  to  the  free  surface.  Figure  14  represents  a  Fictitious  scenario.  The  bubble,  in  this 
case,  would  have  undoubtedly  broken  through  the  free  surface,  venting  its  high-pressure  gases.  It  was 
not  surprising  that  the  calculation  for  this  case  became  unstable  and  inadvertently  terminated. 
Therefore,  understanding  the  limitations  of  the  method  are  important  in  interpreting  the  validity  of  the 
results.  For  a  bubble  at  deeper  depths,  the  bubble  behavior  becomes  less  influenced  by  the  free 
surface.  At  depths  greater  than  100  ft  for  the  300-lb  explosive,  the  bubble  profile  becomes 
indistinguishable  from  the  free  field  profile  for  a  similar  calculatioa 

To  estimate  the  accuracy  of  the  method  for  predicting  reentrant  water  jet  tip  velocities,  a 
calculation  was  compared  with  a  PISCES  code  calculation.  The  PISCES  code  has  been  used  for  a 
number  of  calculations  with  regard  to  underwater  explosion  bubbles,  and  it  is  considered  reasonably 
accurate.  The  comparison  for  peak  reentrant  jet  velocities  was  made  between  the  present  method  and 
PISCES.  If  the  boundary  integral  method  presented  were  valid,  the  two  methods  should  present 
similar  results.  As  it  turned  out,  the  exact  time  and  location  of  the  peak  jet  velocities  varied  slightly 


the  two  methods.  For  the  PISCES  calculation,  the  peak  velocity  occurred  just  behind  the  bubble 
surface  and  was  found  to  be  1,148  ft/s.  A  line  estimating  the  location  of  the  bubble  surface,  with 
regard  to  peak  velocity,  is  shown  in  Figure  15.  Since  an  Euler  method  has  been  used  in  the  PISCES 
calculation,  the  interface  between  different  materials  is  not  defined  exactly.  However,  an  estimated 
profile  is  shown  in  Figure  16.  For  the  incompressible  flow  calculation,  the  maximums  and  minimums 
for  the  potential  must  occur  on  the  boundary  to  satisfy  the  uniqueness  of  Laplace’s  equation.  The 
highest  velocity  recorded  in  the  calculation  occurred  just  as  the  bubble’s  lower  surface  broke  through 
the  bubble’s  upper  surface  (Figure  17).  A  dotted  line  estimating  the  intersection  of  the  two  results  is 
shown  in  Figure  15.  The  results  indicate  a  reasonable  correlation  between  the  two  methods. 

However,  a  more  accurate  estimation  of  the  error  can  only  be  made  through  direct  comparison  with 
experimental  results  which  are  presently  unavailable. 

A  final  calculation  was  performed  for  an  explosion  bubble  near  a  solid  boundary.  The  results 
showing  the  collapse  and  upward  migration  toward  the  solid  boundary  are  summarized  in  Figure  18. 
For  this  calculation,  a  300-lb  charge  of  TNT  was  detonated  12  m  below  the  solid  surface  at  a  depth  of 
12  m.  It  is  interesting  to  note  the  similarities  and  differences  between  the  free  surface  calculation  of 
the  same  weight  and  depth  (Figure  12)  and  this  case.  The  expansion  and  overall  period  of  the  bubble 
growth  are  not  disproportionately  in  disagreement.  However,  as  expected,  the  solid  boundary  increases 
migration  toward  the  boundary  while  having  little  effect  on  the  formation  of  a  bubble  jet.  On  the 
other  hand,  the  free  surface  condition  led  to  decreased  migration  and  a  retardation  of  bubble  jet 
formation.  In  general,  the  free  surface  condition  can  lead  to  bubble  blow-out  (where  the  bubble  vents 
through  the  free  surface),  bubble  downward  migration  at  certain  distances  from  the  free  surface,  and 
the  retardation  of  bubble  jetting.  Therefore,  it  is  clear  that  a  nearby  structure  or  solid  boundary  can 
lead  to  bubble  jetting  into  the  structure  even  when  near  the  free  surface. 

4.  SUMMARY 


The  boundary  integral  method  can  be  used  for  the  prediction  of  bubble  behavior  in  a  variety  of 
scenarios.  While  the  present  study  offers  only  limited  results,  the  method  shows  promise  for  future, 
more  in-depth  studies.  The  method  shows  good  agreement  with  other  analytical  prediction  techniques 
as  well  as  a  variety  of  experimental  data.  This  report  offers  the  only  example  of  such  comparisons 
found  in  the  literature.  The  program  developed  also  offers  a  number  of  improvements  which  include 
efficient  programming  (under  350  lines  of  code)  and  an  accurate  method  for  the  calculation  of  the 


17 


singularity  in  the  problem.  This  approach  allows  the  analytical  study  of  bubble  pulsation  and  jetting 
behavior  in  coordination  with  experimentation.  Future  studies  of  this  nature  should  include  the  three- 
dimensional  aspects  of  bubble  collapse,  thus  bringing  the  method  to  its  full  potential. 

In  summary,  the  methodology  presented  was  simple  to  use,  computationally  efficient,  and 
reasonably  accurate.  The  overall  cost  and  reliability  of  the  method  may  prove  a  valuable  aid  in 
guiding  and  analyzing  experimental  work. 
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Figure  3.  Bubble  Partitioning. 
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Figure  4.  Bubble  Radius  vs.  Time  plot  (.66  lb  TNT  300  ft). 
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Figure  10.  Depth  vs.  Period  Data  (300  lb  TNT). 
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Figure  14.  Bubble  Profiles  (300  lb  TNT  at  6  ml. 
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Figure  16.  PISCES  Calculation  for  Bubble  Profile  0,200  lb  TNT  400-ft  Depth). 
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Figure  18.  Solid  Boundary  Calculation  (300  lb  TNT  12  m). 
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